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Abstract. This paper proposes an algorithmic framework for solving parametric optimization 
problems which we call adjoint-based predictor-corrector sequential convex programming. After 
presenting the algorithm, we prove a contraction estimate that guarantees the tracking performance 
of the algorithm. Two variants of this algorithm are investigated. The first one can be used to 
solve nonlinear programming problems while the second variant is aimed to treat online parametric 
nonlinear programming problems. The local convergence of these variants is proved. An application 
to a large-scale benchmark problem that originates from nonlinear model predictive control of a 
hydro power plant is implemented to examine the performance of the algorithms. 
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1. Introduction. In this paper, we consider a parametric nonconvex optimiza- 
tion problem of the form: 



p(0 



min f{x) 

s.t. g(x)+M£ 
x € ft, 



0. 



where / : R™ — > R is convex, g : R™ — > R m is nonlinear, 51 C R" is a nonempty, 
closed convex set, and the parameter £ belongs to a given subset V C W . Matrix 
M £ M. mxp plays the role of embedding the parameter £ into the equality constraints in 
a linear way. Throughout this paper, / and g are assumed to be differentiable on their 



domain. Problem P(£) includes many (parametric) nonlinear programming problems 



such as standard nonlinear programs, nonlinear second order cone programs, nonlinear 
semidefinite programs [301 133 HS| ■ The theory of parametric optimization has been 
extensively studied in many research papers and monographs, see, e.g. [71 1231 [38] . 

This paper deals with the efficient calculation of approximate solutions to a se- 
quence of problems of the form P(£)[ where the parameter £ is slowly varying. In 
other words, for a sequence {£,k}k>o such that ||£fc+i — £fe|| is small, we want to solve 
the problems P(£fc) in an efficient way without requiring more accuracy than needed 
in the result. 



In practice, sequences of problems of the form P(£) arise in the framework of real 



time optimization, moving horizon estimation, online data assimilation as well as in 
nonlinear model predictive control (NMPC). A practical obstacle in these applications 
is the time limitation imposed on solving the underlying optimization problem for each 
value of the parameter. Instead of solving completely a nonlinear program at each 
sample time [3J SI [27] > several online algorithms approximately solve the underlining 
nonlinear optimization problem by performing the first iteration of exact Newton, 
sequential quadratic programming (SQP), Gauss-Newton or interior point methods 
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[TBI [571 Wf\ . In [ini [37J the authors only consider the algorithms in the framework 
of SQP method. This approach has been proved to be efficient in practice and is 
widely used in many applications [13] . Recently, Zavala and Anitescu |47j proposed 
an inexact Newton-type method for solving online optimization problems based on 
the framework of generalized equations [7j [38] . 

Other related work considers practical problems which possess more general con- 
vexity structure such as second order cone and semidefinite cone constraints, nons- 
mooth convexity 20, 43]. In these applications, standard optimization methods may 
not perform satisfactorily. Many algorithms for nonlinear second order cone and 
nonlinear semidefinite programming have recently been proposed and found many 
applications in robust optimal control, experimental design, and topology optimiza- 
tion, see, e.g. [3] [201 US 123 113] • These approaches can be considered as generalization 
of the SQP method. 

1.1. Contribution. The contribution of this paper is twofold. We start our 
paper by proposing a generic framework for the adjoint-based predictor- corrector se- 
quential convex programming (APCSCP) for parametric optimization and prove a 
main result of the stability of tracking error for this algorithm (Theorem 13. 4p . In 
the second part the theory is specialized to the non-parametric case where a single 
optimization problem is solved. The local convergence of these variants is also proved. 
Finally, we present a numerical application to large scale nonlinear model predictive 
control of a hydro power plant with 259 state variables and 10 controls. The perfor- 
mance of our algorithms is compared with a standard real-time Gauss- Newton method 
and a conventional model predictive control (MPC) approach. 

APCSCP is based on three main ideas: sequential convex programming, predictor- 
corrector path-following and adjoint-based optimization. We briefly explain these 
methods in the following. 

1.2. Sequential convex programming. The sequential convex programming 
(SCP) method is a local nonconvex optimization technique. SCP solves a sequence 
of convex approximations of the original problem by convexifying only the noncon- 
vex parts and preserving the structures that can efficiently be exploited by convex 
optimization techniques [9j [34] [35] . Note that this method is different from SQP 
methods where quadratic programs are used as approximations of the problem. This 
approach is useful when the problem possesses general convex structures such as 
conic constraints, a cost function depending on matrix variables or convex constraints 
resulting from a low level problem in multi-level settings [21 HH S3] • Due to the com- 
plexity of these structures, standard optimization techniques such as SQP and Gauss- 
Ncwton-type methods may not be convenient to apply. In the context of nonlinear 
conic programming, SCP approaches have been proposed under the names sequential 
semidefinite programming (SSDP) or SQP-type methods [ITJ [20j [2TJ [30j [3lJ [43] . It 
has been shown in |17j that the superlinear convergence is lost if the linear semidef- 
inite programming subproblems in the SSDP algorithm are convexificd. In [33j the 
authors considered a nonlinear program in the framework of a composite minimization 
problem, where the inner function is linearized to obtain a convex subproblem which 
is made strongly convex by adding a quadratic proximal term. 

In this paper, following the work in (20j [22] [44] [46], we apply the SCP approach 



to solve problem P(£) The nonconvex constraint g(x) + M£ = is linearized at each 
iteration to obtain a convex approximation. The resulting subproblems can be solved 
by exploiting convex optimization techniques. 
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1.3. Predictor-corrector path-following methods. In order to illustrate the 
idea of the predictor-corrector path- following method [T^l EZ] , we consider the case 
£1 = R n . The KKT system of problem P(£) can be written as F(z;^) — 0, where 
z = (x,y) is its primal-dual variable. The solution z*(£) that satisfies the KKT 
condition for a given £ is in general a smooth map. By applying the implicit function 
theorem, the derivative of z*(-) is expressed as 



(0 



dF 

dz 



dF 



(**(0;0- 



In the parametric optimization context, we might have solved a problem with param- 
eter £ with solution z = 2*(£) and want to solve the next problem for a new parameter 
£. The tangential predictor z for this new solution z*(£) is given by 



z = z*(o + ^(0(i-i) = z*(i) 



dF 

dz 



_1 OF 
~dj 



Note the similarity with one step of a Newton method. In fact, a combination of the 
tangential predictor and the corrector due to a Newton method proves to be useful in 
the case that z was not the exact solution of F(z; £) = 0, but only an approximation. 
In this case, linearization at (z, £) yields a formula that one step of a predictor- corrector 
path-following method needs to satisfy: 



(1.1) 



dF 



dF 



02 



Written explicitly, it delivers the solution guess I for the next parameter £ as 



dF 



(*;£) 



9^ 



<9F 



=Az predictot =Az co „ octor 

Note that when the parameter enters linearly into F, we can write 
^(z;m-0=F(z;0-F(z;0. 



Thus, equation (jl.ip is reduced to 



(1.2) 



F(z;|) + ^(z)(z-z) 



= 0. 



It follows that the predictor-corrector step can be easily obtained by just applying 
one standard Newton step to the new problem P(£) initialized at the past solution 
guess z, if we employed the parameter embedding in the problem formulation |13j . 

Based on the above analysis, the predictor-corrector path-following method only 
performs the first iteration of the exact Newton method for each new problem. In 
this paper, by applying the generalized equation framework [351 132] > we generalize 
this idea to the case where more general convex constraints are considered. When the 
parameter does not enter linearly into the problem, we can always reformulate this 



problem as P(£) by using slack variables. In this case, the derivatives with respect to 
these slack variables contain the information of the predictor term. Finally, we notice 
that the real-time iteration scheme proposed in |16j can be considered as a variant of 
the above predictor-corrector method in the SQP context. 
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1.4. Adjoint-based method. From a practical point of view, most of the time 
spent on solving optimization problems resulting from simulation-based methods is 
needed to evaluate the functions and their derivatives [BJ. Adjoint-based methods 
rely on the observation that it is not necessary to use exact Jacobian matrices of 
the constraints. Moreover, in some applications, the time needed to evaluate all the 
derivatives of the functions exceeds the time available to compute the solution of the 
optimization problem. The adjoint-based Newton-type methods in [18j (26J El] can 
work with an inexact Jacobian matrix and only require an exact evaluation of the 
Lagrange gradient using adjoint derivatives to form the approximate optimization 
subproblems in the algorithm. This technique still allows to converge to the exact 
solutions but can save valuable time in the online performance of the algorithm. 

1.5. A tutorial example. The idea of the APCSCP method is illustrated in 
the following simple example. 

Example. 1.0. (Tutorial example) Let us consider a simple nonconvex parametric op- 
timization problem: 



(1.3) mi 

where £ £ V 
show that x 



{-X! | x\ +2x2+2-4^ = 0, 



1 < 0, x > 0, x e R 2 } , 



K : £ > 1.2} is a parameter. After few calculations, we can 
£ — (2y £ — \/J, 2\/£ — 1) T is a stationary point of problem (|1.3[) which is 
also the uniquely global optimum. It is clear that problem (( 1 . 3|) satisfies the strong 
second order sufficient condition (SSOSC) at xt. 



Note that the constraint x 
second order cone constraint 
M := -4 and Q. := {x £ R 2 | 



casted into the form of P(£) 



(^i,l) T | 



1 < is convex and it can be written as a 
2 < X2- Let us define g(x) := x\ + 2x-2 + 2, 
j < X2, x > 0}. Then, problem (|1.3p can be 




Fig. 1.1. The trajectory of three methods (k = 0, ■ ■ ■ ,9), (o is and o is x k ). 







l|x K -x*(y|| 

SOC const, viol. 











Exact-SQP 



Projected-SQP 



Fig. 1.2. The tracking error and the cone constraint violation of three methods (k = 0, • • ■ , 9). 



The aim is to approximately solve problem (|1.3[) at each given value of the 
parameter £. Instead of solving the nonlinear optimization problem at each until 
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complete convergence, APCSCP only performs the first step of the SCP algorithm to 
obtain an approximate solution x k at . Notice that the convex subproblem needed 
to be solved at each in the APCSCP method is 

(1.4) min {-x l | 2x^xi + 2x 2 -{x\f + 2 - 4£ = 0, || {x u 1) T || < x 2 , x > 0} . 

We compare this method with other known real-time iteration algorithms. The first 
one is the real-time iteration with an exact SQP method and the second algorithm 
is the real-time iteration with an SQP method using a projected Hessian [T6j[29]. In 
the second algorithm, the Hessian matrix of the Lagrange function is projected onto 
the cone of symmetric positive semidefinite matrices to obtain a convex quadratic 
programming subproblem. 

Figures [Ol and HTSI illustrate the performance of three methods when ^ = 1.2 + 
fcA£fc for k = 0, . . . , 9 and = 0.25. The initial point x° of three methods is chosen 
at the true solution of P(£o)- We can see that the performance of the exact SQP 
and the SQP using projected Hessian is quite similar. However, the second order 
cone constraint ||(xi,l) T || 2 < x 2 is violated in both methods. The SCP method 
preserves the feasibility and better follows the exact solution trajectory. Note that 
the subproblem in the exact SQP method is a nonconvex quadratic program, a convex 
QP in the projected SQP case and a second order cone constrained program (11.41) in 
the SCP method. O 

1.6. Notation. Throughout this paper, we use the notation V/ for the gra- 
dient vector of a scalar function /, g' for the Jacobian matrix of a vector val- 
ued function g and S n (resp., 5" and 5+ + ) for the set ofnxn real symmetric 
(resp., positive semidefinite and positive definite) matrices. The notation ||-|| stands 
for the Euclidean norm. The ball B(x,r) of radius r centered at x is defined as 
B(x, r) := {y e M" | \\y — x\\ < r} and B(x, r) is its closure. 

The rest of this paper is organized as follows. Section [2] presents a generic frame- 
work of the adjoint-based predictor- corrector SCP algorithm (APCSCP). Section |3] 
proves the local contraction estimate for APCSCP and the stability of the approxima- 
tion error. Section U considers an adjoint-based SCP algorithm for solving nonlinear 
programming problems as a special case. The last section presents computational 
results for an application of the proposed algorithms in nonlinear model predictive 
control (NMPC) of a hydro power plant. 

2. An adjoint-based predictor-corrector SCP algorithm. In this section, 
we present a generic algorithmic framework for solving the parametric optimization 



problem P(£) Traditionally, at each sample of parameter £, a nonlinear program 
P(£fc) is solved to get a completely converged solution l(£fc). Exploiting the real-time 
iteration idea (T21 [TC] , in our algorithm below, only one convex subproblem is solved 
to get an approximated solution z k at to 2 (£*;)• 

Suppose that z k := (x k ,y k ) £ fix R m is a given KKT point of P(^) (more details 
can be found in the next section), is a given m x n matrix and £ S 1 ?. We 
consider the following parametric optimization subproblem: 

!min {f{x) + (m k ) T (x - x k ) + \{x - x k ) T H k {x - x k )) 
X s.t. A k {x-x k )+g(x k ) + M£ = 0, 
xett, 

where m k := m(z k ,Ak) — (g'(x k ) — Ak) T y k ■ Matrix Ak is an approximation to 
g'(x h ) at x k , Hk is a regularization or an approximation to V 2 £(z k ), where C is the 



G 
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Lagrange function of P(£) to be denned in Section [3] Vector m k can be considered as 
a correction term of the inconsistency between Ak and g'(x k ). Vector y k is referred to 
as the Lagrange multiplier. S ince / and Q are convex and Hk is symmetric positive 



Here, 



Ak and Hk are 



semidefinite, the subproblem P(z k , Ak,-tlk',k,)\ is convex, 
considered as parameters. 

Remark 1. Note that computing the term g'(x k ) T y k of the correction vector m k 
does not require the whole Jacobian matrix g'(x k ), which is usually time consuming to 
evaluate. This adjoint directional derivative can be cheaply evaluated by using adjoint 
methods {2$ . 

The adjoint-based predictor-corrector SCP algorithmic framework is described as 
follows. 



Algorithm 1. (Adjoint-based predictor-corrector SCP algorithm (APCSCP)). 

Initialization. For a given parameter £o G V, solve approximately (off-line) P(£o) 
to get an approximate KKT point z° := (x°,y°). Compute g(x°), find a matrix Aq 

which approximates g'(x°) and Hq € <S" . Then, compute vector m° := (g'(x°) — Aq) T y°. 
Set k := 0. 

Iteration k (k = 0, 1, . . . ) For a given (z k , Ak, Hk), perform the three steps below: 
Step 1. Get a new parameter value £/c+i G T '■ 

Step 2. Solve the convex subproblem P(z k , Ak, Hk; Cfe+i) to obtain a solution 
x k+1 and the corresponding multiplier y k+1 . 

Step 3. Evaluate g(x k+1 ), update (or recompute) matrices Ak+i and Hk+i G S". 
Compute vector m k+1 := g' (x k+1 ) T y k+1 — A^ +1 y k+1 . Set k := k + 1 and go back to 
Step 1. 



The core step of Algorithm [T] is to solve the convex subproblem P(z k , Ak, Hk; £) 
at each iteration. To reduce the computational time, we can either implement an 
optimization method which exploits the structure of the problem or rely on several 
efficient software tools that are available for convex optimization 9, 35, 36]. In this 
paper, we are most interested in the case where one evaluation of g' is very expensive. 
A possibly simple choice of Hk is Hk = for all k > 0. 

The initial point z° is obtained by solving off-line P(£o)- However, as we will show 
later [Corollary 13.5] , if we choose z° close to the set of KKT points Z*(£o) of P(£o) 
(not necessarily an exact solution) then the new KKT point z 1 of P(z°, Aq, Hq; £ ) is 
still close to Z*(£i) of P(£i) provided that ||£i — £ || is sufficiently small. Hence, in 
practice, we only need to solve approximately problem P(£o) to get a starting point 
z°. 

In the NMPC framework, the parameter £ usually coincides with the initial state 
of the dynamic system at the current time of the moving horizon. If matrix Ak = 
g'(x k ), the exact Jacobian matrix of g at x k and Hk = 0, then this algorithm collapses 
to the real-time SCP method (RTSCP) considered in [35]. 

3. Contraction estimate. In this section, we will show that under certain as- 
sumptions, the sequence {z k }k>o generated by Algorithm Q] remains close to the se- 
quence of the true KKT points {zk}k>o of problem P(£fc)- Without loss of generality, 
we assume that the objective function / is linear, i.e. f(x) = c T x, where c £ K" is 



given. Indeed, since / is convex, by using a slack variable s, we can reformulate P(£) 
as a nonlinear program min^s) {s | g(x) + Af£ = 0, x £ il, f(x) < s}. 
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3.1. KKT condition as a generalized equation. Let us first define the La- 
grange function of problem P(£) as 

C(x, y; := c T x + (g(x) + Mtffy, 

where y is the Lagrange multiplier associated with the constraint g{x) + M£ = 0. 
Since the constraint x G ft is convex and implicitly represented, we will consider it 
separately. The KKT condition for P(£) is now written as 



G c + g'(x) T y+Afn{x), 
= g(x)+MZ, 



(3.1) 

where Afn (x) is the normal cone of ft at x defined as 
(3.2) Wn(z) := 



{u G K" | u T (x - v) > 0, v G ft} , if x G ft 
0, otherwise. 



Note that the first line of (|3.1[) implicitly includes the constraint x G ft. 

A pair (x(£), satisfying (|3.1[) is called a KKT point of P(£) and is called 
a stationary point of P(£) with the corresponding multiplier £/(£). Let us denote by 
Z*(£) and A*(£) the set of KKT points and the set of stationary points of P(£), 
respectively. In the sequel, we use the letter z for the pair of (x, y), i.e. z := (x T , y T ) T ■ 

Throughout this paper, we require the following assumptions which are standard 
in optimization. 

Al. The function g is twice differ entiable on their domain. 

A2. For a given £o G V , 'problem P(£o) has at least one KKT point z , i.e. 
Z*(£o)^0- 

Let us define 

(3.3) + 



9{x) 

and K := ft x R m . Then, the KKT condition (13.11) can be expressed in terms of a 
parametric generalized equation as follows: 

(3.4) 0€F(z) + CZ+ff K (z), 

where C := Generalized equations are an essential tool to study many prob- 

lems in nonlinear analysis, perturbation analysis, variational calculations as well as 
optimization [511521155]. 

Suppose that, for some £/. G V, the set of KKT points Z*(£k) of P(£fc) is nonempty. 
For any fixed z k G Z*(^), we define the following set-valued mapping: 

(3.5) L(z;z k ,£ k ) :=F(z k )+F'{z k )(z-z k ) + C(, k +M K {z). 

We also define the inverse mapping L~ l : R n+m — » K ra+m of L(-; z k , as follows: 

(3.6) L- 1 (6;z k ,Z k ):={zeR n+m : (5 G L(z; z k , £/,.)} . 



Now, we consider the KKT condition of the subproblem P(z k , A k , FL k \ if) For 
given neighborhoods B(z k ,r z ) of z k and S(^,r^) of and z fc G B(z k ,r z ), (k+i G 
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B{^ kl r^) and given matrices A k and H k G 5™, let us consider the convex subprob- 
lem P(z k , A k , H k ; ^ k+1 ) with respect to the parameter (z k , A k , H k , £fc+i). The KKT 
condition of this problem is expressed as follows. 



(3.7) 



e c + m(z k , A k ) + H k {x - x k ) + A^y + Mn{x), 
= g(x k ) + A k (x - x k ) + M^ k+1 , 



where Afn(x) is defined by (|3.2[) . Suppose that the Slater constraint qualification 
holds for the subproblem P(z fc , A kl H k \ Cfc+i), i.e.: 

ri(fi) n {x G ffi" | g(x k ) + A k (x - x k ) + M^+i = 0} ^ 0, 

where ri(O) is the relative interior of ft. Then by convexity of fl, a point z k+1 := 
(x k+1 , y k+1 ) is a KKT point of P(z k 7 A k ,H k \ (, k+ i) if and only if x k+1 is a solution to 
P(z fc , Afc, H k ; ^fc+i) associated with the multiplier y k+1 . 

Since g is twice differentiable by Assumption A[T] and / is linear, for a given 
z = (x, y), we have 

n i 

(3.8) V 2 x C(z) = J2yiV 2 9i(x), 

i=l 

the Hessian matrix of the Lagrange function C, where V 2 <?i(-) is the Hessian matrix 
of <j>j (i = 1, . . . , m). Let us define the following matrix: 



H k Al 
A k 



(3.9) F' k 

where H k € 5". The KKT condition (|3.7I) can be written as a parametric linear 
generalized equation: 

(3.10) e F(z k ) + F' k {z~ z k ) + Cti k+1 +N K {z), 

where z k , Ff, and £, k +i are considered as parameters. Note that if A k = g'(x k ) and 
H k = V^,C(z k ) then (|3.10|) is the linearization of the nonlinear generalized equation 
(|3.4|) at (z k ,^ k+ i) with respect to z. 

Remark 2. Note that (|3.10[) is a generalization of (|1.2p , where the approximate 
Jacobian F' k is used instead of the exact one. Therefore, p.lOp can be viewed as one 
iteration of the inexact predictor-corrector path-following method for solving (|3 .4[) . 



3.2. The strong regularity concept. We recall the following definition of the 
strong regularity concept. This definition can be considered as the strong regularity 
of the generalized equation (|3.4[) in the context of nonlinear optimization, see [38j . 

Definition 3.1. Let ^ k e V such that the set of KKT points Z*{£ k ) ofP(^ k ) is 
nonempty. Let z £ Z*(^ k ) be a given KKT point ofP(£, k ). Problem P(£, k ) is said to 
be strongly regular at z k if there exist neighborhoods B(0, fs) of the origin and B(z k , f z ) 
of z k such that the mapping z^(S) := B(z k ,f z ) D L^ 1 (S;z k ,£ ik ) is single-valued and 
Lipschitz continuous in B(0,f$) with a Lipschitz constant < 7 < +oo, i.e. 

(3.11) \\4(s)-4n<-y¥-s'\\, w,<*'eB(o,F 5 ). 

Note that the constants 7, f z and fs in Definition 13.11 are global and do not depend 
on the index k. 
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From the definition of L 1 where strong regularity holds, there exists a unique 
z* k {8) such that 8 G F(z k ) + F'(z k )(z*(S) - z k ) + C£ k + N K (z*(5)). Therefore, 

z* k {5) = (F'^+ATk)- 1 (F'(z k )z k - F(z k ) - C& + 5) 
= J k (F'(z k )z k -F(z k )~C£ k +5), 

where J k := (F'(z k ) +ATk)~ 1 - The strong regularity of P(£) at z k is equivalent to the 
single-valuedness and the Lipschitz continuity of J k around v k := F'(z k )z k — F(z k ) — 

The strong regularity concept is widely used in variational analysis, perturbation 
analysis as well as in optimization 8, 32, 39,. In view of optimization, strong regu- 
larity implies the strong second order sufficient optimality condition (SSOSC) if the 
linear independence constraint qualification (LICQ) holds 38 . If the convex set Q is 
polyhedral and the LICQ holds, then strong regularity is equivalent to SSOSC [T9"] . 
In order to interpret the strong regularity condition of P(£ fc ) at z k G Z*{^ k ) in terms 
of perturbed optimization, we consider the following optimization problem 

min (c - 5 c ) T x + ^(x — x k ) T 'V 2 ,£(x k , y k )(x - x k ) 



s.t. g(x k ) + g'{x k )(x-x k ) + M£ k = S g 



(3.12) 

Here, S = (§ c , S g ) G B(0, f$) is a perturbation. Problem P(£fc) is strongly regular at z k 
if and only if (|3.12l) has a unique KKT point zt(S) in B(z k , f z ) and z k (-) is Lipschitz 
continuous in B(0,f~s) with a Lipschitz constant 7. 

Example. 3.1. Let us recall example (jl.3p in Section 11.11 The optimal multipliers 
associated with two constraints x\ +2x2+2— 4£ = and x\—x\ J r \ < are yl = (2y/£— 

1)[8VC 2 - > and y* 2 = [8\/£ 2 - CV?]" 1 > 0, respectively. Since the last 

inequality constraint is active while x > is inactive, we can easily compute the critical 
cone as C{x\,y*) = {(4,0) G K 2 \ x^di = 0}. The Hessian matrix V 2 x C(x\,y*) = 

2(2/1+1/2) Q £ ^he Lagrange function C is positive definite in C(a;|,y*). Hence, 

the second order sufficient optimality condition for (|1.5p is satisfied. Moreover, y\ > 
which says that the strict complementarity condition holds. Therefore, problem (| 1 . 5[) 
satisfies the the strong second order sufficient condition. On the other hand, it is easy 
to check that the LICQ condition holds for (ll.5[) at x^. By applying [35J Theorem 
4.1], we can conclude that (ll.3|) is strongly regular at (x^,y*). O 

The following lemma shows the nonemptiness of Z* (£) in the neighborhood of £ k . 

Lemma 3.2. Suppose that Assumption A[T1 is satisfied and Z*(£ k ) is nonempty 
for a given G V . Suppose further that problem P(£fc) is strongly regular at z k for 
a given z k G Z*(^ k ). Then there exist neighborhoods B{^ k ,r^) of £ k and B(z k ,r z ) of 
z k such that Z*(£k+i) * s nonempty for all £k+i G B((, k ,r^) and Z*(£ k +i) n B(z k ,r z ) 
contains only one point z k+1 . Moreover, there exists a constant < a < +00 such 
that: 

(3.13) \\z k+1 -z k \\<v\\fai-S k \\. 



Proof. Since the KKT condition of P(£fe) is equivalent to the generalized equation 
(j3.4l) with £ = By applying [38, Theorem 2.1] we conclude that there exist neigh- 
borhoods E(£fc,T"£) of £ k and B(z k ,r z ) of z k such that Z*(£ k+ i) is nonempty for all 
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£,k+i G B(£k, r ?) an d Z*(£k+i) H B(z k , r z ) contains only one point z k+1 . On the other 
hand, since \\F{z k ) + - F{z k ) - C£, k+1 \\ = ||M(& - £ fe+1 )|| < ||M|| - &||, 
by using the formula [351 2.4], we obtain the estimate (|3.13l) . □ 

3.3. A contraction estimate for APCSCP using an inexact Jacobian 
matrix. In order to prove a contraction estimate for APCSCP, throughout this sec- 
tion, we make the following assumptions. 

A3. For a given z k G Z* k > 0, the following conditions are satisfied. 
a) There exists a constant < n < %z such that: 



(3.14) 



F'(z k )-F[ 



where F k is defined by (|3.9I) . 
b) The Jacobian mapping F'(-) is Lipschitz continuous on B(z k ,r z ) around z k , 
i.e. there exists a constant < uj < +oo such that: 



(3.15) 



\F'(z)- F'(z k )\\ <u\\z-z k \\, Vz €B(z k ,r z ). 



Note that Assumption A[3] is commonly used in the theory of Newton-type and 
Gauss- Newton methods [HJ [TS] , where the residual term is required to be sufficiently 
small in a neighborhood of the local solution. From the definition of FL we have 



F'{z k )-F' k 



VlC{z k ) 
g'{x k ) - 



O 



A T k 



Hence, 



F'(z k ) 



pi 



depends on the norms of V^£(z fe ) — H k and g'{x k ) — A k . These 
quantities are the error of the approximations H k and A k to the Hessian matrix 
V 2 x C(z k ) and the Jacobian matrix g'(x k ), respectively. On the one hand, Assumption 
A[3]i) requires the positive definiteness of H k to be an approximation of V 2 ,/^ (which 
is not necessarily positive definite). On the other hand, it requires that matrix A k is a 
sufficiently good approximation to the Jacobian matrix g' in the neighborhood of the 
stationary point x k . Note that the matrix H k in the Newton-type method proposed 
in [7] is not necessarily positive definite. 

Now, let us define the following mapping: 



(3.16) 



J k := {F' k +N K y 



wherc F^ is defined by (|3.9[) . The lemma below shows that J k is single- valued and 
Lipschitz continuous in a neighbourhood of v k := F' k z k — F(z k ) — C^ k . 

Lemma 3.3. Suppose that Assumptions A[TJ A[2] and A[3^) are satisfied. Then 
there exist neighborhoods B(^ k , r%) and B(z k , r z ) such that if we take any z k 6 B(z k , r z ) 
and € B(^ k ,r^) then the mapping J k defined by (|3.16l) is single-valued in a 

neighbourhood B(v k ,r v ), where v k := F' k z k — F{z k ) — C£ k . Moreover, the following 
inequality holds: 



(3.17) 

where (3 := — 



\\J k (v) - J k (v')\\ <@\\v-v'\\, Vv,v' eB(v k ,r v ), 

> is a Lipschitz constant. 

Proof. Let us fix a neighbourhood B(v k , r v ) of v k . Suppose for contradiction that 
Jk is not single- valued in B(v k , r v ), then for a given v the set Jk{v) contains at least 
two points z and z' such that \\z — z'\\ ^ 0. We have 

(3.18) v G F' k z+N K {z) and v G F' k z' + Af K (z r ). 
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Let 
(3.19) 



and 



S:=v- [F' k z k - F(z k ) - C£ k ] + [F'(z ) - F' k ]{z - z k ) 7 



S':=v- [F' k z k - F(z k ) - C&] + [F'(z k ) F' k ](z' - z k ). 
Then (|3.18[) can be written as 



(3.20) 



and 



6 G F(z k ) + F'(z k )(z - z k ) + d k +Af K (z), 



6' G F(z k ) + F'(z k )(z' - z k ) + C^ k +M K {z'). 
Since v in the neighbourhood B(v k , r v ) of v k := F k z k — F(z k ) — C£ fc , we have 

\\5\\< \\v-v k \\+ [F'{z k )-Fl](z-z k ) 



<r v + 



F\z k )-F' k 



\z — z 



< r v + n\\z — z . 

From this inequality, we see that we can shrink B(z k ,r z ) and B(v k ,r v ) sufficiently 
small (if necessary) such that \\S\\ < fg. Hence, S G 23(0, fa). Similarly, 5' G 23(0, fg). 
Now, using the strong regularity assumption of P(£fc) at z k , it follows from (I3.20P 

that 



(3.21) \\ Z - Z '\\< 7 \\ S -S'\\. 
However, using (|3.19p . we have 

\\5-6'\\= \[F'(z k )~F k }(z-z') 
< F'{z k )-F' k \\z-z' 



era 
< k||z-z'||. 

Plugging this inequality into (|3.21l) and then using the condition jk < i < 1, we get 

\\z-z'\\ <\\z-z'\\, 

which contradicts to z ^ z'. Hence, J k is single- valued. 

Finally, we prove the Lipschitz continuity of Jk- Let z = Jk(y) and z' — J k (v'), 
where v, v' G B(v k , r v ). Similar to (13.201) . these expressions can be written equivalently 
to 



(3.22) 
where 
(3.23) 



and 



6 G F(z k ) + F'(z k )(z - z k ) + CZ k +Mk(z), 



and 



6' G F(z k ) + F'(z k )(z' - ~z k ) + C£ k +M K (z'), 

5:=v- [F' k z k F(z k ) C£„] + [F'{z k ) - F' k ]{z - z k ), 
6' := v 1 - [F' k z k - F(z k ) - C&] + [F'(z k ) - F' k ]{z' - z k ). 
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By using again the strong regularity assumption, it follows from (|3.22p and (|3.23[) 
that 



u*-/ii<7iis-<ni 



< 



7 || v — v || +7 



[F'{z k )-F' k ](z-z') 



< 7 ||i> — v || + 7K ||2 — z || . 
Since jft < ^ < 1, rearranging the last inequality we get 



\z-z'\\ < 



1 — 7« 



\\v-v'\\, 



which shows that satisfies (|3. 17|) with a constant (3 



> 0. □ 



Let us recall that if z h+1 is a KKT of the convex subproblem P(z k , Ak, iffc;6+i) 



then 



e F' k (z k+1 - z k ) + F{z k ) + C6+1 +H K {z k+1 ). 



According to Lemma f3 .31 if z k S B(z k ,r z ) then problem P(z k , A)., H^; £) is uniquely 
solvable. We can write its KKT condition equivalently as 



(3.24) 



z k+1 = J k (F' k z k - F(z k ) - C6+1 



Since z k+1 is the solution of (|4.2[) at 6+1 ) we have = F(z 
where u k+1 e MK{z k+1 ). Moreover, since 



-,k+i 



fc+i) . 



k+l 



C6+i + u 
we can write 



(3.25) 



i k+1 = Jk (Kz k+1 - F(z k+1 ) - C&+: 



The main result of this section is stated in the following theorem. 

Theorem 3.4. Suppose that Assumptions A[T1-A[2l are satisfied for some 6 € P '■ 
Then, for k > and z k G Z*(6) 7 */P(6) * s strongly regular at z k then there exist 
neighborhoods B(z k ,r z ) and B(^k,T^) such that: 

a) The set of KKT points Z*(6+i) of P(6+i) is nonempty for any 6+1 £ 
S(6,r e ). 

b) If, in addition, Assumption A[3a) is satisfied then subproblem P(z k , A/,, H k \ 6+i) 
is uniquely solvable in the neighborhood B(z k ,r z ). 

c) Moreover, if, in addition, Assumption A[3b) is satisfied then the sequence 
{z k }k>o generated by Algorithm^ where 6+i £ £?(£&, r^), guarantees 



y k+i 



(3.26) 



~ k z k \ 



J fe+1 ||< (a + Cl ||z fc -z fc ||) \\z k 

+ (C2 + C 3 ||6+1 - 611) 116+1 - 611 , 



where < a < 1, < Q < +oo, i = 1, . . . , 3 and ci > are given constants 



and z k+1 e Z*(6+i). 



Proof. We prove the theorem by induction. For k = 0, we have Z*(£o) is nonempty 
by Assumption A[2j Now, we assume Z*(6) is nonempty for some k > 0. We will 
prove that Z*(6+i) is nonempty for some 6+i £ ^(Ck,r^), a neighborhood of 6- 

Indeed, since Z*(6) is nonempty for some 6 £ Pi we t & ke an arbitrary z k £ 
Z*(6) such that P(6) is strong regular at z k . Now, by applying Lemma 13.21 to 
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problem P(6)j then we conclude that there exist neighborhoods B(z k ,r z ) of z k and 
B{£k,rt) of £ fc such that Z*(£ fc+1 ) is nonempty for any £ fc+1 £ B(^ k ,r e ). 

Next, if, in addition, Assumption A[3^) holds then the conclusions of Lemma |3~3" 
hold. By induction, we conclude that the convex subproblem P(z k , A k , £&) is uniquely 
solvable in B(z k , r z ) for any £ k+ i S B(£ k , r^). 

Finally, we prove inequality (13.261) . From (|3.24l) . (13.251) and the mean- value the- 
orem and Assumption A[3}d), we have 



I fe+ i_„ fe+ i||l3p 



J k ( {Flz k - F(z k ) - d k+l ) - z 



(3.27) 



< P 



= (3 



J k (^z fe - F{z k ) - Ca+i) - J k (F' k z k+1 - F(z k+1 ) - a 

F' k {z k - z k+1 ) - F(z k ) + F(z k+1 ) 
F' k {z k -z k )-F{z k )+F{z k )\ + \F(z k+1 )~F(z k )~F^(z k+1 -z k ) 

[F k - F'(z k )](z k — z k ) — f [F\z k +t(z k -z k ))-F'(z k )](z k -z k )dt 



\Fl-F\z k )}{z k+1 -z k ) - / [F'{z k +t(z k+1 -z k ))-F'(z k )}(z k+1 -z k )dt 



< 



\z k+1 ~z k 



By substituting (|3~T3)l into (|3~27l) we obtain 

"ll<K-+|ll^ 



\z k+1 -z k+1 



\z k -z k \ 



-P ( «cr+ — 116+1 



-611 116+1 -61 



and 



If we define a := 0k = ^ < 1 due to A[3ji), ci := > 0, c 2 := > 

c 3 := 2(1 -7k) — as f° ur given constants then the last inequality is indeed (|3.26l) . □ 
The following corollary shows the stability of the approximate sequence {z k } k >o 

generated by Algorithm Q] 

Corollary 3.5. Under the assumptions of Theorem \3.4\ there exists a positive 

number < r z < f z :— (1 — ajc^ 1 such that if the initial point z° in Algorithm^]] is 

chosen such that \\z° — z°|| < r z , where z° 6 Z*(£o) then, for any k> 0, we have 

(3.28) \\z k+1 - z k+1 \\ <r z , 

provided that ||6+i _ 611 ^ r £; where z k+1 6 Z*((, k +i) and < < f £ with 



(2C3)" 1 + 4c 3 r z (l - a - Cl r z ) - c 2 

c 2 _1 r z (l - a- c x r z ) 



ifc 3 > 0, 
«/c 3 =0. 



Consequently, the error sequence {e k } k >o, where e k :— ||z fc — z* , between the exact 
KKT point z k and the approximate KKT point z k ofP(^ k ) is bounded. 
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Proof. Since < a < 1, we have r z := (1 — a)c 1 1 > 0. Let us choose r z such that 
< r z < f z . If z° e B(z°,r z ), i.e. ||z° - z°\\ < r z , then it follows from (j3"T2!))) that 

jjz 1 ~z 1 \\<{a + Cl r z )r z + (c 2 + c 3 - Coll) 116 — Co 1 1 - 

In order to ensure [jz 1 — z 1 ]] < r z , we need (c 2 + C3 — fo||) ||£i — foil < P : = 
(1 — a — c ir z )r z . Si nce < r z < f z , p > 0. The last condition leads to ||£i — fo|| < 
(2c 3 )- 1 ( v /cf+ 4c 3 p-c 2 ) if c 3 > and ||gi - go|| < - a - ar x ) if c 3 = 0. By 

induction, we conclude that inequality (I3.28[) holds for all k > 0. □ 

The conclusion of Corollarv l3.5l is illustrated in Figure |3~TI where the approximate 
sequence {z k }k>o computed by Algorithm Q] remains close to the sequence of the true 
KKT points {z k }k>o if the starting point z° is sufficiently close to zq. 




£0 Si in 5fc+i 

Fig. 3.1. The approximate sequence {z fc }fc>o along the trajectory z(-) of the KKT points. 



3.4. A contraction estimate for APCSCP using an exact Jacobian ma- 
trix. If Ak = g' (x k ) then the correction vector m k = and the convex subproblem 
P(z k , Ak, Hk] £) collapses to the following one: 



{c T X+±(x-x k ) T H k (x-x k )\ 



mm \c x 

s.t. g(x k )+g'(x k )(x-x k ) + M£ 
x e Q. 



0, 



Note that problem P(x k ,Hk]£,) does not depend on the multiplier y k if we choose 
.Hfc independent ly of . We refer to a variant of Algori thm [1] where we use the con- 



vex subproblem |P(z fc , Hk; £) instead of P(z fc , Afc, ijfc; £)| as a predictor- corrector SCP 
algorithm (PCSCP) for solving a sequence of the optimization problems {P(£fc)}fc>o- 
Instead of Assumption A[3^,) in the previous section, we make the following as- 
sumption. 



A3'. There exists a constant < k < ^- such that 

— 27 



\V 2 x C{z k ) - H k \\ <K,Vk> 0. 



(3.29) 



where V^£(z) defined by 

Assumption A[33requires that the approximation Hk to the Hessian matrix W^C(z k ) 
of the Lagrange function C at z k is sufficiently close. Note that matrix Hk in the 
framework of the SSDP method in [TT] is not necessarily positive definite. 
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Example. 3. 5. Let us continue analyzing example (| 1 . 3[) . The Hessian matrix of the 
Lagrange function C associated with the equality constraint x\ + 2x2 + 2 — 4£ = 



is V|£(z|, yi ) 



2vl 




, where y* is the multiplier associated with the equality 



constraint at x|. Let us choose a positive semidefmite matrix Hk := ['"q 1 o]' wnere 



hn > 0, then 



VlC(x*yl)-H k 



we can choose hu > such that 



\Vt — /in | - Since j/^ > 0, for an arbitrary k > 0, 
— J/J] < «• Consequently, the condition (I3.29[) 
is satisfied. In the example (II .3|) of Subsection II .51 we choose hu = 0. O 
The following theorem shows the same conclusions as in Theorem 13.41 and Corol- 
lary [3?5] for the predictor- corrector SCP algorithm. 

Theorem 3.6. Suppose that Assumptions A[T1-A[2l are satisfied for some £o G V . 
Then, for k > and z k £ Z*(£k), «/P(£/c) is strongly regular at z k then there exist 
neighborhoods B(z k ,r z ) and fi(£/-,rf) smc/i </mt: 

a) The set of KKT points Z*{^k+i) of P(£fc+i) is nonempty for any £fc+i G 
B(&,r € ). 

&,) //, in addition, Assumption A|3H is satisfied then subproblem P(x k , H k ; £k+i) 
is uniquely solvable in the neighborhood B(z k ,r z ). 

c) Moreover, if, in addition, Assumption A[3b) then the sequence {z k } k >o gen- 
erated by the PCSCP, where £k+i € B(^k,r^), guarantees the following in- 
equality: 



zk+1 I 



(3.30) 



where < a < 1, < c\ < 
and z k+1 G 



< (a + ci ||;T - z K ||) - 

+ (p2 + C8||&+l-&||)||ffc+l-&||, 

+oo, i = 1, • • • , 3 and £2 > are given constants 



d) If the initial point z° in the PCSCP is chosen such that \\z° — z°|| < r z , where 
z° G Z*(£ ) and < f z < r z :— cj" 1 (l - a), then: 



(3.31) ||z fc+1 -z K+1 || <f„ 

provided that ||£fc+i — £fc|| < ^ wii/i < f £ < r^, 



=fe+l| 



r 5 



(2c 3 ) 



V^l 



y/% + 4c 3 r x (l 
-a — c\f z ) 



cif z ) - c 2 



ifh > 0, 

«/C3 =0. 



Consequently, the error sequence {\\z k — z k \\}k>o between the exact KKT 
point z k and the approximation KKT point z o/P(£fc) is still bounded. 
Proof. The statement a) of Theorem 13.61 follows from Theorem l3.4l We prove b). 
Since Ak = g'(x k ), the matrix F' k defined in (|3.9[) becomes 

h . = \ H k g'(x k j 
k ' [g'{x k ) J ' 

Moreover, since g is twice differentiable due to Assumption A[Tl g' is Lipschitz contin- 
uous with a Lipschitz constant L g > in B(x k ,r z ). Therefore, by Assumption AI3 7 ! 
we have 

2 r w2 />f-k\ _//=fc\T -.//„ife\T1 2 

F'(z fc ) 



V*£(* fc ) g , {x k ) T -g'(x k ) T - 
g'(x k )-g'(x k ) 



(3.32) 



< ||V2£(z fc )-ijJ| 2 + 2|| ff '( a; fc )- ff '(a; fc )| 
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Since £7 < |, we can shrink B(z k , r z ) sufficiently small such that 



If we define hi := \/k 2 + 2L 2 r 2 > then the last inequality and p. 321) imply 



(3.33) 



F'(z k 



where £17 < tj- Similar to the proof of Lemma 13.31 we can show that the mapping 
Jfc := (FL+Nk) is single- valued and Lipschitz continuous with a Lipschitz constant 
P := 7(1 — 7Ki) _1 > in B(z k , r z ). Consequently, the convex problem P(x k , Hk',£k+i) 
is uniquely solvable in B{z k ,r z ) for all £k+i € B(£,k,r^), which proves b). 



With the same argument as the proof of Theorem 13. 4| we can also prove the 
following estimate 



v fc+i 



z k \\ < (a k + ci \\z k -z fe ||) \\z k ~z k \\ + (c 2 + c 3 ||&+i ~&l 



where a :— 7«i(l — 7%) g [0,1), Si := 7^(2 — 27K1) > 0, £2 := 7Ki(t(1 — 
l7«i) _1 > and c 3 := jLua 2 (2 — 2"fki)~ 1 > 0. The remaining statements of Theorem 
13.61 are proved similarly to the proofs of Theorem 13.41 and Corollary 13.51 □ 

Remark on updating matrices Ak and Hk ■ In the adjoint-based predictor-corrector 
SCP algorithm, an approximate matrix Ak of g'(x k ) and a vector m k — (g'{x k ) — 
Ak) T y k are required at each iteration such that they maintain Assumption A[3j Sup- 
pose that the initial approximation Aq is known. For given z k and Ak, k > 0, we 
need to compute Ak+i and m k+1 in an efficient way. If \\Ak — g'(x k+1 )\\ is still small 
then we can even use the same matrix Ak for the next iteration, i.e. Ak+i = Ak due 
to Assumption A[3] (see Section [5]). Otherwise, matrix Ak+i can be constructed in 
different ways, e.g. by using low-rank updates or by a low accuracy computation. As 
by an inexactness computation, we can either use the two sided rank-1 updates (TR1) 
PS [26] or the Broyden formulas [S] . However, it is important to note that the use of 
the low-rank update for matrix Ak might destroy possible sparsity structure of matrix 
Ak . Then high-rank updates might be an option [51 125) . 

In Algorithm [T] we can set matrix Hk = for all k > 0. However, this matrix 
can be updated at each iteration by using BFGS-type formulas or the projection of 
V 2 x C{z k ) onto SI. 

4. Applications in nonlinear programming. If the set of parameters E col- 
lapses to one point, i.e. E := {£} then, without loss of generality, we assume that 
£ = and problem P(£) is reduced to a nonlinear programming problem of the form: 



(P) 



min fix 
s.t. g(x) = 0, 



T 

c x 



where c, g and f2 are as in P(£) In this section we develop local optimization algo- 
rithms for solving (|P|) . 

The KKT condition for problem ([P]) is expressed as: 



(4.1) 



€ c + g'{x) T y + Nn{x), 
= g(x), 
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where A/q(x) defined by (|3.2I) . A pair z := (x T ,y T ) T satisfying (|4.ip is called a 
KKT point, x is called a stationary point and y is the corresponding multiplier of 
(|P]) , respectively. We denote by Z* the set of the KKT points and by S* the set of 
stationary points of (|Pl) . 

Now, with the mapping F defined as (13.31) . the KKT condition (|4.1[) can be 
reformulated as a generalized equation: 

(4.2) OeF(z)+M K (z), 

where K = £1 X M m as before and JVk (z) is the normal cone of if at z. 
The subproblcm P(z fc , Afc, H^; £) in Algorithm [1] is reduced to 

!min c T x + {jmP) T (x — x J ) + ^(x — x^) T Hj[x — x J ) 
"s.t. 5 (x J ') + A,-(x- x^') = 0, 
x e n. 

Here, we use the index j in the algorithms for the nonparametric problems (see below) 
to distinguish from the index k in the parametric cases. 

In order to adapt to the theory in the previous sections, we only consider the 
full-step algorithm for solving which is called full-step adjoint-based sequential 
convex programming is described as follows. 



Algorithm 2. (Full-step adjoint-based SCP algorithm (FASCP)) 



Initialization. Find an initial guess x° € fi and y° G 

mated to g'(x°) and H G 5". Set m° := ( 5 '(x° 



A ) T y° and j 



a matrix Aq approxi- 
0. 



Iteration j. For a given (z 3 , Aj, Hj), perform the following steps: 



Step 1. Solve the convex subproblem P(z 3 , Aj, Hj) to obtain a solution x{ +1 and 



the corresponding multiplier y 



Step 2. If 



J+i 



< e, for a given tolerance e > 0, then: terminate. Other- 



wise, compute the search direction Ax J 



S'tep 5. Update x 



i+i — 



Evaluate the function value g(x J+ ), update 



(or recompute) matrices j4j+i and -ffj+i G 5™ (if necessary) and the correction vector 



Increase j by 1 and go back to Step 1. 



The following corollary shows that the full-step adjoint-based SCP algorithm generates 
an iterative sequence that converges linearly to a KKT point of (|Pj) . 

Corollary 4.1. Let Z* ^ and z* G Z*. Suppose that Assumption AfT] 
ZioMs onrf that problem (jPj is strongly regular at z* (in the sense of Definition \3. j]) . 
Suppose further that Assumption A[3K) *s satisfied in B(z*,r z ). Then there exists a 
neighborhood B(z*,r z ) of z* such that, in this neighborhood, the convex subproblem 
P(x : ' , Aj, Hj) has a unique KKT point z^ +1 for any z^ G B(z*,r z ). Moreover, if, in 
addition, Assumption A[3b) holds then the sequence {z J }j>o generated by Algorithm 
[H starting from z a £ B(z* , r z ) satisfies 

(4.3) \\z j+1 - z*\\ <(a + c 1 ||z J - z*\\) \\z j - z*\\ , Vj > 0, 

where < a < 1 and < C\ < +00 are given constants. Consequently, this sequence 
converges linearly to z* , the unique KKT point of (|Pj in B(z*,r z ). 

Proof. The estimate (|4.3I) follows directly from Theorem l3.4l bv taking ^ = for 
all k. The remaining statement is a consequence of (14.31) . □ 
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If Aj = g'{x^) then the convex subproblem P(z- J , Aj, Hj) in Algorithm[2]is reduced 

to: 



mm 



n {c T x+\{x-x j ) T Hj(x-x j )} 

P^'.-Hj) \ s.t. g(x:>)+g'(x j )(x-x j )=0, 

x e ft. 

The local convergence of the full-step SCP algorithm considered in |44j follows from 
Theorem 13.61 as a consequence, which is restated in the following corollary. 

Corollary 4.2. Suppose that Assumption A^holds and problem (|P]) is strongly 
regular at a KKT point z* e Z* (in the sense of Definition \3.i\) . Suppose further that 
Assumptions ^[3^1 and fi^p) are satisfied. Then there exists a neighborhood B(z*,r z ) 
of z* such that the full-step SCP algorithm starting from x° with (x , y ) G B(z* , r z ) 
generates a sequence {z J }j>o satisfying: 

< (ti + Ci^J - z*\\) Wz 1 - z*\\ , 

where a G [0, 1) and Cj G [0, +oo) are constants and z J+1 is a unique KKT point of 
the subproblem P(x J , Hj). As a consequence, the sequence {z 3 } converges linearly to 
z* , the unique KKT point of JFJ) in B(z* ,r z ). 

Finally, it is necessary to remark that if ft is a polyhedral convex set in R n , i.e. 
f2 is the intersection of finitely many closed half spaces of K™ , then problem fP]) also 
covers the standard nonlinear programming problem. It was proved in [19 that if f2 
is polyhedral convex and the constraint qualification (LICQ) holds then the strong 
regularity concept coincides with the strong second order sufficient condition (SSOSC) 
for fP]). In this case, by an appropriate choice of Hk, the SCP algorithm collapses 
to the constrained Gauss-Newton method which has been widely used in numerical 
solution of optimal control problems, see, e.g. [6]. 

5. Numerical Results. In this section we implement the algorithms proposed 
in the previous sections to solve the model predictive control problem of a hydro power 
plant. 

5.1. Dynamic model. We consider a hydro power plant composed of several 
subsystems connected together. The system includes six dams with turbines Di (i = 
1, . . . , 6) located along a river and three lakes L\, L2 and L3 as visualized in Fig. 15.11 
U\ is a duct connecting lakes L\ and L2. T\ and Ti are ducts equipped with turbines 
and Ci and C2 are ducts equipped with turbines and pumps. The flows through the 
turbines and pumps are the controlled variables. The complete model with all the 
parameters can be found in [5D] . 

The dynamics of the lakes is given by 

,,,, dh(t) _ g in (t) - g ou t ft) 

{5A > ~dT~ S ' 

where h(t) is the water level and S is the surface area of the lakes; qi n and q ou t are the 
input and output flows, respectively. The dynamics of the reaches Ri {i = 1, . . . , 6) is 
described by the one-dimensional Saint- Venant partial differential equation: 



(5.2) 



9q(t,y) , ds(t,y) _ n 
dy " r dt u ' 



gdt \s(t,y)J + 2gdy \ S Ht,y)J + dy +1 f{ l >y) M^J ~ U. 
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Fig. 5.1. Overview of the hydro power plant. 



Here, y is the spatial variable along the flow direction of the river, q is the river flow (or 
discharge), s is the wetted surface, h is the water level with respect to the river bed, 
g is the gravitation acceleration, If is the friction slope and la is the river bed slope. 
The partial differential equation (|5.2j) can be discretized by applying the method of 
lines in order to obtain a system of ordinary differential equations. Stacking all the 
equations together, we represent the dynamics of the system by 

(5.3) w(t) = f(w,u), 

where the state vector w <E R n ™ includes all the flows and the water levels and 
u e R"" u represents the input vector. The dynamic system consists of n w — 259 states 
and n u — 10 controls. The control inputs are the flows going in the turbines, the 
ducts and the reaches. 

5.2. Nonlinear MPC formulation. We are interested in the following NMPC 
setting: 



min J(w(-),u(-)) 

w.u 

(54) s.t. w = f(w,u), w(t) = w (t), 

u(t) £ U, w{t) eW, TE[t,t + T] 
w{t + T) e R T , 

where the objective function J(w(-),u(-)) is given by 

rt+T 

J(w(-),u(-)) := / [(w(t) - w s ) t P(w(t) - w s ) + (u(r) - u s ) t Q(u(t) - u s )] dr 
(5 5) 

+ (w(t + T) - w s ) T S(x(t + T) ~ w s ). 

Here P, Q and S are given symmetric positive definite weighting matrices, and (w s , u s ) 
is a steady state of the dynamics (|5.3I) . The control variables are bounded by lower 
and upper bounds, while some state variables are also bounded and the others are 
unconstrained. Consequently, W and U are boxes in and 1"» , respectively, but 
W is not necessarily bounded. The terminal region Rt is a control-invariant ellipsoidal 
set centered at w s of radius r > and scaling matrix S 1 , i.e.: 

(5.6) R T := {w e W* w \ (w - w s ) T S(w - w s ) < r} . 
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To compute matrix S and the radius r in (|5.6|) the procedure proposed in |10) can be 
used. In [25J it has been shown that the receding horizon control formulation (|5.4I) 
ensures the stability of the closed-loop system under mild assumptions. Therefore, 
the aim of this example is to track the steady state of the system and to ensure the 
stability of the system by satisfying the terminal constraint along the moving horizon. 
To have a more realistic simulation we added a disturbance to the input flow qi n at 
the beginning of the reach R\ and the tributary flow ^tributary- 
The matrices P and Q have been set to 



P := diag 



Q := diag 



0.01 



(«>«)? + 1 



1 < i < n,, 



{ui + u b )f + 1 



1 < i < n, 



where it; and Ub is the lower and upper bound of the control input u. 

5.3. A short description of the multiple shooting method. We briefly 
describe the multiple shooting formulation [B] which we use to discretize the continu- 
ous time problem (|5.4p . The time horizon [t, t + T] of T — 4 hours is discretized into 
Hp = 16 shooting intervals with Ar = 15 minutes such that tq — t and Tj+i := r.; + Ar 
(i = 0, . . . , H p — 1). The control u(-) is parametrized by using a piecewise constant 
function u(t) = Ui for n < r < Tj + Ar [i = 0, . . . , H p — 1). 

Let us introduce H p + 1 shooting node variables (z = 0, . . . , iip). Then, by inte- 
grating the dynamic system w — f(w, u) in each interval [n, T{ + Ar], the continuous 
dynamic (|5.3[) is transformed into the nonlinear equality constraints of the form: 



(5.7) 



g{x) + M£ := 



so - £, 
w(s 0l u ) - Sl 

W(S H „-l,«H „-l) - S Hp 



= 0. 



Here, vector x combines all the controls and shooting node variables ui and Si as 
x = (sj, iip, . . . , Sft _y, Ujj _i, s^) T , ^ is the initial state wo(t) which is considered 
as a parameter, and w(ui,Wi) is the result of the integration of the dynamics from r; 
to Ti + Ar where we set u(t) — m and w(ji) = Sj. 
The objective function (|5.5I) is approximated by 



H„-l 



f(x):= ^ [(si-w s ) T P(si-w s ) + (ui-u s ) T Q(i 



0] 



(5.8) 



+(sh p - w s ) T S(s Hp - w s ), 

while the constraints are imposed only at t — Ti, the beginning of the intervals, as 

(5.9) G W, m 4 € 17, s Hp G R T , (i = 0,.-.,H p - 1). 

If we define O := U Hp x (VK^p x Rt) C M™ x then is convex. Moreover, the objective 
function (|5.8I) is convex quadratic. Therefore, the resulting optimization problem is 



indeed of the form P(£) Note that Q is not a box but a curved convex set due to Rt- 
The nonlinear program to be solved at every sampling time has 4563 decision 
variables and 4403 equality constraints, which are expensive to evaluate due to the 
ODE integration. 
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5.4. Numerical simulation. Before presenting the simulation results, we give 
some details on the implementation. To evaluate the performance of the methods 
proposed in this paper we implemented the following algorithms: 

• Full-NMPC - the nonlinear program obtained by multiple shooting is solved 
at every sampling time to convergence by several SCP iterations. 

• PCSCP - the implementation of Algorithm[T]using the exact Jacobian matrix 
of g. 

• APCSCP - the implementation of Algorithm Q] with approximated Jacobian 
of g. Matrix Ak is fixed at Ak = g'(x°) for all k > 0, where ir is approximately 
computed off-line by performing the SCP algorithm (Algorithm [2]) to solve 



the nonlinear programming P(£) with £ = £o = wo(t). 
• RTGN - the solution of the nonlinear program is approximated by solving 
a quadratic program obtained by linearizing the dynamics and the terminal 
constraint sh € Rt- The exact Jacobian g'(-) of g is used. This method can 
be referred to as a classical real-time iteration [16] based on the constrained 
Gauss- Newton method [51 . 
To compute the set Rt a mixed Matlab and C++ code has been used. The computed 
value of r is 1.687836, while the matrix S is dense, symmetric and positive definite. 

The quadratic programs (QPs) and the quadratically constrained quadratic pro- 
gramming problems (QCQPs) arising in the algorithms we implemented can be ef- 
ficiently solved by means of interior point or other methods [§J 135) . In our imple- 
mentation, we used the commercial solver CPLEX which can deal with both types of 
problems. 

All the tests have been implemented in CH — h running on a 16 cores workstation 
with 2.7GHz Intel®Xeron CPUs and 12 GB of RAM. We used CasADi, an open 
source C++ package pQ which implements automatic differentiation to calculate the 
derivatives of the functions and offers an interface to CVODES from the Sundials pack- 
age [42j to integrate the ordinary differential equations and compute the sensitivities. 
The integration has been parallelized using openmp. 

In the full-NMPC algorithm we perform at most 5 SCP iterations for each time 
interval. We stopped the SCP algorithm when the relative infinity-norm of the search 
direction as well as of the feasibility gap reached the tolerance e = 10 -3 . To have a 
fair comparison of the different methods, the starting point x° of the PCA, APCA 
and RTGN algorithms has been set to the solution of the first full-NMPC iteration. 

The disturbance on the flows qi n and ^tributary are generated randomly and varying 
from to 30 and to 10, respectively. All the simulations are perturbed at the same 
disturbance scenario. 

We simulated the algorithms for H m = 30 time intervals. The average time 
required by the four methods is summarized in Table 15.11 Here, AvEvalTime is the 



Table 5.1 
The average time of four methods 



Methods 


AvEvalTime [s] 


AvSolTime [s] 


AvAdjDirTime [s] 


Total [s] 


Full-NMPC 
PCSCP 
RTGN 
APCSCP 


220.930 (91.41%) 
70.370 (90.05%) 
70.588 (96.97%) 
0.458 ( 3.28%) 


20.748 (8.58%) 
7.736 (9.90%) 
2.171 (2.98%) 
11.367 (81.34%) 


2.122 (15.18%) 


241.700 
78.142 
72.795 
13.975 



average time in seconds needed to evaluate the function g and its Jacobian; AvSolTime 
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is the average time for solving the QP or QCQP problems; AvAdjTime is the average 
time for evaluating the adjoint direction g'(x k ) T y k in Algorithm[T] Total corresponds 
to the sum of the previous terms and some preparation time. On average, the full- 
NMPC algorithm needed 3.27 iterations to converge to a solution. 

It can be seen from Table lOl that evaluating the function and its Jacobian matrix 
costs 90% — 97% of the total time. On the other hand, solving a QCQP problem is 
almost 3 — 5 times more expensive than solving a QP problem. The computationally 
expensive step at every iteration is the integration of the dynamics and its lineariza- 
tion. The computational time of PCSCP and RTGN is almost similar, while the time 
consumed in APCSCP is about 6 times less than PCSCP. 

The closed-loop control profiles of the simulation are illustrated in Figures 15.21 
and 15.31 Here, the first figure shows the flows in the turbines and the ducts of lakes 
L\ and L2, while the second one plots the flows to be controlled in the reaches Ri 
(i = 1, . . . , 6). We can observe that the control profiles achieved by PCSCP as well 








- - - rtscp 




rtgn 




rtadjscp 


1 - 
1 (__! 


full-scp 







Fig. 5.2. The controller profiles Qt-d 9Ci> Qt 2 and qc±- 

as APCSCP are close to the profiles obtained by Full-NMPC, while the results from 
RTGN oscillate in the first intervals due to the violation of the terminal constraint. 
The terminal constraint in the PCSCP is active in many iterations. 

Figure 15.41 shows the relative tracking error of the solution of the nonlinear pro- 
gramming problem of the PCSCP, APCSCP and RTGN algorithms when compared 
to the full-NMPC one. The error is quite small in PCSCP and APCSCP while it 
is higher in the RTGN algorithm. This happens because the linearization of the 
quadratic constraint can not adequately capture the shape of the terminal constraint 
sn 6 Rt- The performance of APCSCP is nearly as good as PCSCP. This feature 
confirms the statement of Corollarv l3.5l 

6. Conclusions. We have proposed an adjoint-based predictor- corrector SCP 
algorithm and its variants for solving parametric optimization problems as well as 
nonlinear optimization problems. We proved the stability of the tracking error for the 
online SCP algorithms and the local convergence of the SCP algorithms. These meth- 
ods are suitable for nonconvex problems that possess convex substructures which can 
be efficiently handled by using convex optimization techniques [45] . The performance 
of the algorithms is validated by a numerical implementation of an application in non- 
linear model predictive control. The basic assumptions used in our development are 
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Fig. 5.3. The controller profiles of qR x , . . . , qn 6 . 
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Fig. 5.4. T/ie relative errors of PCSCP, APCSCP and RTGN compared to Full-NMPC. 



the strong regularity, Assumption A[3Jd) and Assumption A[3ji) (or Al^]). The strong 
regularity concept introduced by Robinson in [38j and is widely used in optimization 
and nonlinear analysis, Assumption A[3b) (or is needed in any Newton-type 

algorithm. As in SQP methods, these assumptions involve some Lipschitz constants 
that are difficult to determine in practice. 

Our future work is to develop a complete theory for this approach and apply it to 
new problems. For example, in some robust control problem formulations as well as 
robust optimization formulations, where we consider worst-case performance within 
robust counterparts, a nonlinear programming problem with second order cone and 
semidefinite constraints needs to be solved that can profit from the SCP framework. 
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